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I. 



INTRODUCTION 



The current design trend toward higher tip speed, 
reduced engine weight, and higher pressure ratios has made 
the modern high speed fan and compressor more susceptible to 
supersonic unstalled flutter. Unstalled flutter, as 
distinguished from stall flutter, occurs at small angle of 
attack and with attached flows. In addition, it occurs at or 
near the design condition. Flutter itself is a self-excited, 
unsteady phenomenon which can occur when the exciting forces 
just equal the elastic forces of the turbomachinery blades. 
For a review of the problem of fan and compressor flutter, 
refer to Ref. 6. 

The original work in flutter analysis evolved from the 
flutter of aircraft wings, which was observed as early as 
1914. Theodorsen [Ref. 22] developed the most basic analytic 
treatment of unsteady, subsonic aerodynamics for an airfoil 
oscillating in pitch and plunge. This analysis dealt with 
what has become known as classical flutter, and showed that 
flutter was only possible when the natural frequencies of 
flexure and torsion coalesced as the relative velocity over 
the wing was increased. Each mode of oscillation by itself 
(pitch or plunge) was stable. 

Cascade flutter differs from classical flutter in that 
the inertia and stiffness properties of a turbomachine blade 
differ greatly from those of a wing. As a result, the 
critical flutter speed, the speed at which the onset of 
flutter occurs, is very much greater. Garrick and Rubinow 
[Ref. 9] extended the work of Theodorsen to include 
supersonic flow, and their analysis predicted that single 
degree of freedom torsional flutter was possible, although 
the bending mode was still stable. Lane [Re,:. 14] showed 
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that cascade flutter could be analyzed by considering only a 
single oscillating blade. 

Work on cascades with subsonic leading edge locus is 
fairly recent. The first attempt was made by Gorelov 
[Ref. 10]. Verdon [Ref. 23] obtained some of the first 
numerical results using a finite difference procedure. Brix 
and Platzer [Refs. 3,4] were also able to obtain results 
using a method of characteristics approach. Another 
solution was presented by Nagashima and Whitehead [Ref. 15] 
using dipole distributions. All three solutions were in 
close agreement. Sisto and Ni [Ref. 20] developed a 'time 
marching" technique applicable to an infinite cascade. Their 
results appear to be in good agreement with the finite 
cascade results. 

An analytic development, valid for low frequency blade 
motion, has been developed by Kurosaka [Ref. 11] and is 
currently being extended to the higher frequency range 
[Ref. 12]. These results are applied to an infinite cascade. 
Very recently, another analytic treatment has been presented 
by Verdon [Ref. 25] (see also Verdon and McCune [Ref. 24]) 
in which an infinite series representation is used to 
express the influence of neighboring blades and wakes on the 
reference passage. However, series convergence was not 
always obtained, and some doubt exists as to its uniform 
validity [Ref. 13]. A new analysis currently being 
developed by Verdon, applicable to the infinite cascade, 
does appear to yield uniformly good results [Ref. 26]. 
Actual wind tunnel tests of cascade flutter have been made 
by Snyder and Commerford [Ref. 19] and by Fleeter [Ref. 7]. 

In the current study, a linearized method of 
characteristics procedure, based on a developement by Teipel 
[Ref. 21] for unsteady aerodynamics of a flat plate- 
airfoil, is used to examine the pressure distributions and 
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flutter boundaries for cascade blades having subsonic 
leading edge locus. It is an extension of previous work done 
by Chalkley [Ref. 5] and Brix [Ref. 2]. 

In addition, the method of characteristics was applied 
to the problem of supersonic flow through oscillating 
cylindrical shells to determine the pressure distribution 
and associated generalized forces acting on the cylinder 
walls. 

Although it is unreasonable to suppose that the 
two-dimensional cascade investigation can accurately model 
the complex mixed transonic-supersonic flow of an actual 
present day turbomachine, it is hoped that this analysis 
will provide a better understanding of the problems 
involved, and serve as a starting point for the flutter 
analysis during the design stage of modern high speed fans 
and compressors. 
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II 



THEORY 



A. PROBLEM FORMULATION 

Consider a cascade composed of a semi-infinite array of 
thin, two-dimensional flat plate airfoils immersed in the 
supersonic flow of an inviscid, non-conducting, ideal gas 
having constant specific heats. Additionally, the flow 
field is assumed to be irrotational and isentropic. 



The only restriction on the cascade is that the locus of 
blade leading edges must be subsonic; i.e.. 



Y 




/ 



/ 




X 



Figure 1. Cascade geometry 




(II- 1) 
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where /3 is the stagger angle and a is the Mach angle 

1 

( sin- 1 - 1 — ). As shown in figure 1, the cascade is aligned m 
M 

the X-Y coordinate system so that the leading-edge of the 
first blade is at the origin and the blade itself lies along 
the X-axis. 

Perturbations in the flow field are caused by blade 
motion in pitch and plunge, assumed to be of small amplitude 
with simple harmonic motion. For a cascade with subsonic 
leading-edge locus, the blades are swept aft of the initial 
Mach wave. Thus, disturbances are allowed to propagate 
throughout the domain of influence (the region downstream of 
the Mach cone emanating from the leading edge of the first 
blade) , and are not restricted to the region between blade 
passages. 



B. GOVERNING EQUATIONS 



Subject to the initial assumptions stated above, the 
governing eguations for this boundary value problem are the 
continuity eguation 



Dp du dv 

— + P { — + — ) 

DT dX dY 



0 



the eguation of motion ( Euler equations ) 

Du + I dP _ Q 
DT p dx 

— ♦ M . 0 

DT p dY 

the condition of irrotationality 



(II- 2) 



(II-3a) 



(II-3b) 
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0 



(II- 4} 



du dv 
dY dX 



and the energy eguation 



Ds 



DT 



0 



(II- 5) 



D 



denotes the material, or substantive, derivative of a 



fluid particle which is a function of position and time. It 

is composed of the local derivative (rate of change with 
respect to time) and the convective derivative (rate of 
change with respect to position) and is written as 



where V is the velocity vector. The local sonic velocity is 
given by. 



DT 



dp s dp 



(II- 6) 



and 



P 



(II- 7) 



constant 



Taking the total differential of Eg. (7) 



i y t (r’) 

(j) d P - yp(j) dp = o 



(II- 8) 



Rearranging 




(II- 9) 



or 
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(11-10) 



c 2 




P 



C. BOUNDARY CONDITIONS 

The flow boundary condition (flow tangency) for a blade 
having an equation of the form 

F (X,Y,T) = 0 (11-11) 

is developed as follows. At each point of a solid-fluid 
surface, at every instant, the normal component of the 
relative velocity between the solid and the fluid must 
vanish. In ether words, the total time rate of change in 
F (X , Y , T) following a surface particle of the solid is zero. 
Mathematically, this statement is expressed as 

DF 

— =0 (II-12) 

DT 

For the cascade, the equations for the upper and lower 
surfaces of a blade are written 

F (X, Y ,T) = Y - Y (X , T) = 0 (H-13) 

u u 

and 



F^(X,Y,T) = Y - Y^(X,T) 



Applying Eq. (12) to Eqs. (13) 



DF U dl 0 _dF u dFy 

— ■ = — u + u— + v — = 0 
DT dT ax a Y 



= 0 

and (14) 



(11-14) 



(11-15) 



or 



v 



dY u _dY u 
u — - 0 

aT ax 



(11-16) 
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since 



1. Likewise 



d£u 

dY 



DF. dY, dY, 

— L = v L - u— = 0 

DT dT dX 


(H-17) 



The other boundary condition that must be enforced is the 
Sommerfeld radiation condition, i.e., disturbances must 
propagate away from their sources. 

D. LINEARIZATION 

The normal velocity perturbation is written 



dY u - dY u 
v (X, Y (X,T) ) = — + U -— 1 1 

u d T dX 


(11-18) 


dY, dy L 

v(X, Il (X.T,) = U-L 


(11-19) 



From the assumption of thin blades and small amplitude 
oscillations, all flow quantities may be considered to be 



small quantities which are linearly superimposed 
freestream quantities. Thus 


onto the 


u = u + u 

0 


(11-20) 


c = c + c 

0 


(11-21) 


V = V 

while the density and pressure perturbations are 


(11-22) 


Ap = p - p Q 


(11-23) 
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A p = p - p 



(11-24) 



Also, because of the assumption of thin blades, Y and Y 

u L 

are small quantities with respect to axial distances along 

the blade. Thus, substituting Eqs. (20) and (22) into Egs. 
(18) and (19) and neglecting the effect of higher order 
terms, the normal velocity perturbation equation becomes 



v (X, Y (X,T) ) 
u 



a Yu dY u 
— + u — 
d T <>dX 



(11-25) 



v(X,Y l (X,T)) 



dY, a Y, 

— - + u — - 
dT odx 



(11-26) 



To transfer the boundary condition to the axis, Eqs. 
(25) and (26) are expanded in a Taylor series about Y = 0. 



v(X,Y (X,T) ) 
u 



dv (X , Y=0' s ‘ , T) 

- v (X,Y = 0 + ,T) + Y 

u dY 

1 d 2 v (X,Y = 0 + ,T) 

-Y 2 + e © «. 

2 u (dY) 2 



(11-27) 



v (X, Y^ (X,T) ) = V(X,Y = 0-,T) + Y^ 

1 d 2 v (X,Y = 0~,T) 

— Y 2 - 

2 L (dY) 2 

Neglecting higher order terms, 

v (X, Y (X, T) ) = v (X,Y = 0 + ,T) 



dv (X,Y=0-,T) 



dY 



(11-28) 



(11-29) 



v(X,Y (X,T)) = v(X,Y=0~,T) 
and thus 



(11-30) 
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V (X, 0+,T) 



(11-31) 



<rt_u + _ai M 

ai °ax 



v (X, 0~,T) 



ax, ay, 

— + u — - 

ai °ax 



(11-32) 



th 

For cascaded blades, it is assumed that the a blade 

st 

leads the motion of the (n-1) blade by an amount 8, the 

interblade phase angle. Expressing this and the assumption 
of sinusoidal motion, the equation for the upper and lower 
blade surfaces becomes 



Y (X, T) = Y (X) e 



i[wT + (n- 1 ) 8 ] 



(11-33) 



Y^ (X,T) = Y^ (X) e 



i[wT ♦ (n-1 ) 8 J 



(11-34) 



which reduces the flow tangency condition, Eqs. (31) and 
(32), to 



aY i £ (n- 1)8} iwT 

v(X,0+,T) = [i w Y (X) + u ]e e 

u °ax 



(11-35) 



ay i { (n- 1 ) 8} icuT 

v (X, 0~,T) = [iwY (X) + u ]e e 

L L ' <>dX 



(11-36) 



Now let Y (X, T) = Y (X,T) = Y (X) e 

u L 



Then the flow 



tangency condition for pitch is. 



Y (X) = - a (X - x c) 

pitch o o 



(11-37) 
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(11-38) 



itoT 

a = a e 
o 

while for the plunge. 



Y (X) 

plunge 



where. 



iwT 

h = h e 
o 



h £ 
o 



(11-39) 



(11-40) 



Deflections are defined as positive downward, and angle of 

twist positive leading edge up. h and a are dimensionless 

o o 

complex amplitudes, with x q the twist axis location (elastic 

axis) in percent chord c. Figure 2 shows the positive 
direction of all quantities used in defining the flow 
tangency condition. 



y 




Figure 2. Positive directions of Lift 
and Moment resulting from deflections 
in pitch (a) and/or plunge (h) . 



To non-dimensionalize the flow tangency equations, 
distances are scaled to the blade chord and time to blade 
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chord divided by freestream velocity u q . This leads to the 
introduction of the reduced frequency parameter 



(OC 

u G 


(11-41) 



which is twice the reduced frequency as defined by Garrick 
and Rubinow [Ref. 9] who scaled all quantities to blade 
semi-chord. 

Non-dimensionalizing Egs. (35) and (36) and using the 
definition of reduced frequency above, the flow tangency 
conditions become 
pitch: 



v (x,0±, t) 

U 0 


ikt 

= - a Q {cos[ (n-1)8] - k(x-X Q )sin[ (n-1)8]}e 




ikt 



- ia Q £ksin[ (n- 1) 8 ] + k (x-x q ) cos[ (n- 1) S ]} e 



plunge: 


(11-42) 


v (x,0±,t) 
U o 


ikt 

= h Q (ksin[ (n- 1) 8] - ikcos[ (n- 1) 8 ]} e 




(11-43) 



To linearize the expression for the pressure, Eqs. (10) , 
(20) and (24) are used. Expand Eg. (10) in terms of 
perturbation quantities. Through Eq. (48) , second order 
terms are neglected whenever they occur. 



c2 + 2cc 
0 0 


(P 0 + Ap) 

■ 7 (1 ♦*£, 
P 
o 



A p 

As -y- < 1 , the denominator can be expanded in a geometric 
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series of the form 



1 

1 + € 

to yield. 



€ + € 2 “ € 3 + 



(11-45) 



C 2 + 2cc = y(p + Ap) (1 - 
0 0 1 ' r o r# p 

Expanding, 



(11-46) 



c 2 Ap 

c 2 + 2cc = (1 e -~ -) A p 

0 0 y A p 

However, as Ap and Ap are both small quantities 



Sa. A’ 1 = J- 

Ap dp c 2 



(11-47) 



(11-48) 



Expanding Eg. (47) , substituting into 

rearranging terms, 

Ap 2 2 yc 

= cc [ 1 + ] 

P 0 r- 1 0 ( 7 - 1 )c 0 



(46) and 



(11-49) 



2 yC 

But < 1. Using Eg. 

(r-i)c 0 

brackets yields the linearized 
difference 



(43) to expand the term in 
form of the pressure 



P 



7-1 



pc (c - c ) 



(11-50) 
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III. METHOD OF CHARACTER 1ST ICS 



A. CHARACTERISTIC DIRECTIONS 

To introduce the method of characteristics, it is 
important to recall that for supersonic flow there are 
distinct regions of influence and dependence. Hence, curves 
must exist along which the values of the flow variables and 
their derivatives are known, yet the flow field in the 
immediate neighborhood of these curves is indeterminate. 
Therefore, to obtain the directions in the x-y plane of 
these curves, which are called characteristic directions, or 
more simply, characteristics, one must prescribe that the 
derivatives across these curves are indeterminate. 

Mathematically, this means that the dependent variables 
must be continuous, but their first, or higher order, 
derivatives may be discontinuous. Physically, the 
characteristics act as information carriers, in that ail 
disturbances propagate along them. In particular, they 
represent the paths of Mach waves, and are therefore 
referred to as Mach lines. 

To facilitate the analysis, it is desirable to introduce 
new coordinates £ and ^ . The governing equations are 
written in terms of the two arbitrary, but intersecting 
curves 

£ = £ (x,y) = constant (Ill-la) 

and 

V = V (x,y) = constant (Ill-lb) 
However, the equations are more suitably transformed to the 
arbitrary coordinates if the dependent variables are u, v 
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and <3. Here, air is assumed to be a perfect gas with 



constant specific heats c (constant pressure) 

P 

(constant volume). Rewriting Eg. (11-10) as 


and c 

V 


c 2 = yR6 


(III- 2) 


R = c - c 
P V 


(III- 3) 


where R is the gas constant for air, and 


9 is 



temperature in degrees Rankine. Taking the total 
differential of Eg. (2) and dividing by c 2 , 



dc d9 

c 9 

Taking the total differential of Eg. (11-10) , 


(III- 4) 


1 d p 

2cdc = v-dp - p — 
1 P P z 


(III- 5) 


or, dividing by c 2 , 

dc dp dp 

c P p 


(III- 6) 



Now from the First Law of Thermodynamics and the fact that 
the flow is homentropic. 



p 

9ds = du dp 

int p z 


(III- 7) 



For an ideal gas, du = c de and Eg. (7) may be written 
int v 

as 



dp 

9ds = c d9 - R9 — 

V p 


(III- 8) 
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Dividing by c 0 



v 



1 dG d p 

— ds = —• - (y-1)—- (III- 9) 

c y e p 

and using Egs. (4) and (6) yields 



1 dc 1 

— ds = 2- (y-1) -dp (III-10) 

c v c p 



But because the thermodynamic properties of a continuous 
flow field are to be considered as continuous, single- valued 
scalar point functions. Eg. (10) must hold for any change of 



the fluid properties s, c and p. Hence, it may be written in 
terms of material derivatives. Then, subject to Eg. (H-5) , 
Eg. (10) becomes 



1 Ds 1 Dc 

c DT c DT 

v 



(/-I) 



1 Dp 
pDT 



0 



Then, using Eg. (11-10), Eg. (11) reduces to 

Dp 2 _Dc 

DT /-I DT 



(III-11) 



(III-12) 



Now the perturbation eguations are introduced. Writing Eg. 
(II-9) as 



Dp Dp 

= c2— 

DT DT 








(III- 13) 


and substituting 


into Eg 


. (12) allows Eg. 


(II-2) to be 


expressed as 










2 dc 2 

y-1dT y- 


dc 

-u — — + 
1 °dX 


du 

c 

°d X 


dv 

+ c — — = 0 

odr 


(III-14) 


neglecting the 


effect 


of 


higher order terms. This is the 



form of the continuity eguation which is to be used in the 
transformation. 
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To render the first Euler equation, Eq. (H-3) , in a 
form suitable for transformation, use Eq. (11) in 

differential form 



Y _ 1 

2— dc = (r-l)-dp 
c p 

With Eq. (II-9) this yields 


(III-15) 


1 2 

— dp = cdc 

p r- 1 


(III-16) 



Expanding the total derivatives and equating the partial 



derivatives, one obtains 




I^P 2 dc 

/>dX " /-I dX 

Substituting into the first Euler equation 


(III- 17) 


Du du _du du 2 _dc 

= — - + u— + v— = - c ~ — 

DT dT dX dY Y - 1 dX 


(III-18) 



Introducing the small perturbation quantities and neglecting 
second order quantities yields the desired form of the first 
Euler equation. 



du du 2 dc 

- + u — + c — = 0 

dT odx /-I °dx 


(III-19) 


The condition of irrotationality 

d v du 

dY ' dx = ° 


(III-20) 



is used in place of the second Euler equation as one of the 
governing equations to be transformed. Then the continuity 
equation, the first Euler equation and the condition of 
irrotationality form a system of three equations in the 

three unknowns u, v and c. For the assumption of simple 
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harmonic motion, the Teipel amplitude functions [Ref. 21] 
are introduced. 



iwT 

U (X,Y) e 




(III-21a) 



icoT 

V (X , Y) e 



(III-21b) 



icoT 2 1 c - c 

C(X,Y)e = ( — ) (III-21C) 

y-i V c o 

where U, V and C are complex non-dimensional amplitudes. U 
and V are the perturbation velocities and C is the 
perturbation in sonic velocity which will be related to the 
pressure coefficient through Eg. (II--49) . 



Substituting Egs. (21) into Egs. (14), (19) and (20), 
and then non-dimensionalizing as before yields, 

dU , dv dc 

— + /M2-1 + M2— + ikM 2 c = 0 (III-22a) 

ox dy ox 



dO dc 

— + — + ikU = 0 (III-22b) 

dx dx 



dU , d V 

/M2-1 — =o (III-22c) 

dy dx 



For £ = constant or i) = constant, both £ (x,y(x)) and 
yiXtYix)) are implicit functions of x. Hence, 



d£ = £ dx + £ dy = 0 
x y 



and 



(III-23a) 
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(III-23b) 



d?? = 17 dx 






0 



Therefore , 




and the labeling of 
arbitrary. By convention 
normally denoted £ and 



(III- 24) 

the characteristic directions is 
the left-running Mach line is 
the right-running Mach line, ^ . 



To solve for the characteristic directions in the x-y 
plane, Eqs. (22) are rewritten along lines of constant £ 
and rj . 

£ U + /M2-1 £ V ♦ M 2 £ C = - t) U - M z-n C - 

x £ y £ x £ x 17 x 17 

/m 2 -1 -n V - ikM 2 c (III-25a) 

y 77 



£0 + £ C = - <« U -fl C - ikU (III-25b) 

x £ ^x £ 'x 'x 77 



£ U - £ V = - 77 U + /P-i v V (III-25C) 

y x y 77 x 77 



Now Cramer's rule is applied to solve for the 
derivatives of the dependent variables with respect to £ 
along 77 = constant. To satisfy the condition that these 

0 

derivatives be finite requires that each be of the form — , 

0 

i.e., indeterminate. Equating the determinant of the 

t 

coefficient matrix of ( U , C , V } to zero yields, 

£ [ £2(M 2 ~1) - £2 l = 0 (III- 26) 

xx y 
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to which there are three real solutions, each of which 
defines a characteristic direction in the x-»-y plane. These 
are. 



dy 

dx 



e 

x 



0 



(III-27a) 



dy 

dx 



(III-27b) 



The same results are obtained when one solves for the first 
derivatives with respect to 77 along lines of £ = constant. 



Eg. (27a) gives the streamline slope in the flow field, 
while Eg. (27b) gives the slope of the two curves £ = 
constant and 77 = constant. It should be noted that these 
are the same characteristic directions as obtained for 
steady, two-dimensional flow. By convention. 



dy 1 

( — ) = + ; 

dx £ =const /M 2 - 1 



(II I-2Qa) 



dy 1 

dx rj=const /fi 2 - 1 



(III-28b) 



ft 

dx stream 



(III-28C) 



For an arbitrary function f = f (x,y) , the total 

derivative is 

df df 

df = — dx + dy (III- 29) 

dx dy 

For f(x,y) taken along a curve as defined by Egs. (1) , 
however, y is an implicit function of x and Eg. (29) is 
rewritten as 
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(Ill- 30) 



df _ df df dy 

dx dx dydx 

Therefore, along a streamline, 

df _ df 

^dx^ stream dx 

while along the Mach lines. 



df 


df 


1 


df 


dx £ =const 


dx 


+ /m 2 - 1 


dy 


df 


df 


1 


df 


dx rj =const 


dl 


/m 2 -1 


dy 



Then, solving for the partial derivatives in 
ordinary derivatives, 

df df 

T- = (7-) v 
dx dx str 



d£ _ 2 ££ df 

dx 2 dx £ =const dx 77 =const 



df 1 , df df 

— = — /H2-1 £ ( ) _ ( ) ] 

dy 2 dx £ =const dx 77=const 



The compatibility relations are obtained by 
(22) using Egs. (33) . Thus, 

1 dU dU 1 dV dV 

2 dx £ dx 7] 2 dx £ dx rj 

M2 dC dC 

+ “ [ (— ) ♦ (— ) 3 + ikM2 C = 0 

2 dx £ dx 77 



(III- 31) 

(III-32a) 

(III-32b) 
terms of the 

(III-33a) 

(II I-33b) 

(III-33c) 

writing Egs. 

(III-34a) 
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1 dll 
-[ (— ) 

2 dx £ 



dU 1 dC 

W * 1 [ Wc 



dC 

(— ) 3 

dx 7 ] 



ikU 



(III-34b) 



dU dU dV dV 

C (— K - (T-) 3 - C (T-), + (— ) 3 = 

dx £ dx 7 ) dx ( dx 77 



(III-34C) 



Eg. (34b) is multiplied by M 2 and then subtracted from Eg. 
(34a) . The resulting equation is first added to, and then 
subtracted from Eg. (34c) to give two of the compatibility 
equations 

dU dV M 2 

( — ) - ( — ) + ik (U - C) = 0 (III-35a) 

dx £ dx £ M 2 -1 

dU dV M 2 

( ) + ( — ) * ik (U - C) = 0 (III-35b) 

dx 17 dx 77 M 2 -1 

The third compatibility relation is obtained from writing 
Eg. (34b) along a streamline. 



dU dC 

( — ) + ( — ) + ikU = 0 (III-35C) 

dx stream dx stream 



The governing equations have now been reduced to a system of 
three ordinary differential equations which are implicitly a 
function of x, with U , V and C complex amplitudes. Thus, 
Egs. (35) are really six eguations, with three real and 
three imaginary components. 



To write the boundary conditions in terms of the Teipel 
amplitude functions, note that 



v (x, 0 ,t) 



u 

0 



V (x, y) e 



ikt 



(III- 36) 



from the definition of the amplitude function, and therefore 
the flow tangency condition is. 
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V = V (x,y) = 



1 v(x,0,t) 

prf 7" 



(III- 37) 



v (x,0 , t) 

where is given by (11-40) for pitch and (11-41) for 

u o 

plunge. 



B. FINITE DIFFERENCE EQUATIONS 



To write Egs. (III-34) in finite difference form, the 
computational molecule shown in figure 3 is used. 

y 




x 



Figure 3. Computational molecule for a 
general flow field point. 

If F is an arbitrary complex flow quantity, derivatives 
with respect to x along the characteristics, in finite 
difference fcrrn, are written, 

dF - F, , 

(— ) = LL . (III-38a) 

dx str Ax 
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(III-3eb) 



dF 

(— ) c 
dx £ 




( H) = F 22 ~ F 2I 

dx V ^ AX 



The values of the perturbation quantities in 
are averages, so that 



(F) = — ( F ♦ F ) 

str 2 11 22 



(F) = — ( F F ) 

£ 2 12 22 



(F) = — ( F + F } 

7 ] 2 22 21 



1 . G en er al Flow Fie ld Eq uati ons 

With these relations substituted into Eqs. (35 ) , 
of equations becomes, 

U +C -AU = K 
22R 22R I 221 12R 



U +C + A U = K 
221 221 I 22R 121 



U -V - B ( U -C )=K 
22R 22R I 221 221 34R 



U -V + B ( U -C ) = K 

221 221 I 22R 22R 341 



U + V - B ( U -C ) = K 
22R 22R I 221 221 56R 



(III-38C) 

Eqs. (38) 
(III-39a) 
(III-39b) 
(III-39c) 

the system 
(III'40a) 
(III-40b) 
(III-40C) 
(III-40d) 
(III-40e) 
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u 


+ 


V 


+ 


B ( U 




- 


c ) 


= K 




(II 1-4 Of ) 


221 




221 




I 


22 R 


22R 


561 




K 




U 


+ 


C 


♦ 


a a 








(III-40g) 


12R 




1 1R 




1 1 R 




i 


111 








K 


= 


U 


+ 


C 




A U 


1 1 R 






(III-4Gh) 


121 




111 




111 




I 








K 


= 


U 




V 


+ 


B ( 


U 


C 


) 


(III-40i) 


34R 




12R 




12R 




I 


121 


121 




K 


= 


U 




V 




B ( 


U 


C 


) 


(III-40 j) 


341 




121 




121 




I 


12R 


12R 






K 




U 


+ 


V 


+ 


B ( 


U 


C 


) 


(III-40k) 


56R 




2 1 R 




21 R 




I 


211 


211 




K 


= 


U 


+ 


V 




B ( 


U 


C 


) 


(III-401) 


561 




211 




211 




I 


2 1 R 


2 1 R 



and 



1 

A = -k Ax 


(III-4 la) 


I 2 





1 M 2 

B = -k( ) A x ‘ (III-41b) 

I 4 M 2 -1 

As shown by Teipel [Ref. 21], solving these equations gives 



( 1 K 34 R tK 56 R ) ~ B I K 12T ] * 2B l[2 ( K 34I* K 56T } * B r K 1?fJ 



(1-AjBj) 2 + ( 2Bj ) 2 



(I I I- 42a) 



(1 WW- B,^] - 2Bjl(K, 4R , 



(I II- 42b) 



V 

22R 



= — ( K 

2 56R 



34R 



(III-42c) 
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) 



(III-42d) 



V = — ( K - K 
221 2 561 341 



C = K -U + A U 
22R 12R 22R I 221 



(III-42e) 



C = K - U -AU 
221 121 221 I 22R 



(II I- 42f ) 



These are the finite difference equations for a grid point 
located in the general flow field. 

2. U££er Surf ace of a B lade 



For a grid point on top of a blade, the 
computational molecule of figure 3 is modified as shown in 
figure 4. 

y 




Figure 4. Computational molecule for a grid 
point on the upper surface of a blade. 

Perturbation quantities at grid point F are considered 

12 

zero and the normal velocity (pitch or plunge at grid 

point F 2^) is prescribed by the blade movement as given by 
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Eg. (11-37) . Here the applicable equations are 



u 


♦ 


C 


A U 


K 




(III-43a) 


22R 




22R 


I 221 


12E 






u 


+ 


C + 


A U 


K 




(III-43b) 


2 21 




221 


I 22R 


121 






u 




B 0 


+ B C 


= K 


- K 


(III-43C) 


22R 




I 221 


I 221 


56R 


34 R 




u 


+ 


B U 


- B C 


= K 


- K 


(III-43d) 


221 




I 22R 


I 22R 


561 


341 




where K 




and K 


are as 


before i 


r but 




12 




56 








K 


_ 


V 








(III-43e) 


34E 




22R 










K 




V 








(III-43f) 


341 




221 











as given by the flow tangency condition. Blade position is 
always given relative to the leading edge of each blade. 
Solving, 



u = (l ' A r B I l [ K 56R~ K 34R~ 8 I K 12l] >2B l[ K 56I ~ *341* B I K 12 r] 



’-wk 



i k i ;r]~ 






(II I -44a) 



(III -44b) 



C = K - U + A U 
22R 12R 22R I 221 



(III-44C) 



C = K - U - A U 
221 121 221 I 22R 



and V are known. 
22 



(III-44d) 
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3 . Lower Surface of a Blade 



At a grid point on the lower surface of a blade, the 
normal perturbation velocity at grid point is also 

known. The computational molecule used is shown in figure 



y 




x 



Figure 5. Computational molecule used for a 
grid point on the lower surface of a blade. 



For this grid point, the finite difference equations 

become. 



[ 


+ 


C 


- A U 


= 


K 




(III-45a) 


22R 




22R 


I 


221 


12R 




1 


+ 


C 


+ A U 


_ 


K 




(III-45b) 


2 21 




221 


I 


22R 


121 




1 


_ 


B ( U 


C 


) = K 


+ K 


(III-45c) 


22R 




I 


221 


221 


56R 


34R 


1 


+ 


B ( U 


C 


) = K 


+ K 


(III-45d) 


221 




I 


22R 


22R 


561 


341 


K 

12 




is the same 


as Egs . (39a) 


and (39b) , 


but 
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(II I-46a) 



K - U 

56E 12R 

K = U 

561 121 

K = V 

34R 22R 

K = V 

341 221 



V ‘ ♦ B (U - C ) 

12R I 121 121 



V - B (U - C ) 
121 I 12R 12R 



(III-46b) 

(III-46C) 

<III-46d) 



Solving for the perturbation quantities at grid point 22 
yields. 



U 22R 



( 1 - ) [k 56R , K 34r - B 1 K 12I ]t2B I [K 56I ♦ K 34I + 

( 1 - AjBj ) 2 + (2Bj) 2 

(| - A I B 1> f K 56I ^ *341 ^ B I K 12 r1- 2B i['W K 34R- B l K 12l] 

M - AjBj ) 2 ♦ ( 2 B j) 2 



(I I I-47a ) 



(III- 47b) 



C = K 

22R 12R 

C = K 

221 121 

and V are known. 



U + A U 

22R I 221 

U - A U 

221 I 22R 



(III-47C) 

(III-47d) 



4. In itial Left R un ning Mach Line 



To obtain the equations for the initial left running 
Mach line ( £ = constant) , the computational molecule of 
figure 5 is used. However, it is assumed that grid points 

F and F are just in the freestream and the limiting 

11 21 

condition x-*-0 is taken. For x-*-0, A^_ and B^-^0. Thus Egs. 
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(38) reduce to. 



U + V =0 
22 22 



(III-48a) 



U + C =0 
22 22 



(III-48b) 



or. 



U = - V = - C 
Now Eg. (35a) becomes 



(III- 49) 



dU M 2 

( — ) = - ik U 

dx £ M2-1 



(III- 50) 



which can be integrated along £ = constant. Or, 



U = U| expf-ik x] 

I x=0 M2-1 



(III- 51) 



U I is known from Egs. (49) and (37) Thus, 

I x=0 

U = - V (O)cos(k x) - V (O)sin(k x) 

22E 22R M2-1 221 M2-1 



(III-52a) 



M 2 M 2 

U = V (0) sin (k x) - V (0) cos (k x) 

221 22R M 2 - 1 221 M 2 -1 

(III-52b) 

with the appropriate form of Eg. (37) used for V (0). Then, 

22 



V = - U 
22R 22R 



(III-52C) 



V = - U 
221 221 



(III-52d) 



C = - U 
22R 22R 



(III-52e) 



45 



(III-52£) 



C 



U 



221 



221 



Egs. (52) are also used to compute the perturbation 
quantities at the leading edge of the upper surface of the 
first blade. 

5- Init ial Right Ru nn ing Ma ch Line 



along the initial right running Mach line (77 = constant) is 

analagous to the development for the left. Here. F and F 

12 11 

of figure 2 are assumed to be just in the freestream. Eqs. 
(38) reduce to, 

U - V =0 (III-53a) 

22 22 



The developement for the perturbation quantities 




(III-53b) 



or. 



U = V 



c 



(III- 54) 



Then Eq. (35b) becomes. 



dU 




(III- 55) 



Integrating along 77 = constant yields 



U 




(III- 56) 



or. 




M2 

( 0 ) sin (k x) 

M 2 - 1 



(III-57a) 
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U = - V (O)sin(k x) 

221 22R M2-1 



Again, the appropriate form of Eg. 



V (0) . Then, 
22 



V = U 
22R 22R 



M2 

+ V (O)cos(k x) 

221 M2-1 

(II I-57b) 

(37) is used to define 



(III-57c) 



V =0 
221 221 



(III-57d) 



(III-57e) 



221 



221 



Eqs. (57) are used to compute the 

the first grid point on the 

blade. In the computer program, 

grid point F are subscripted 33 
22 

respectively. 

6. First Grid Point on a New 



(III-57f) 

perturbation quantities at 
lower surface of the first 
lower surface points at 
or 55 for pitch and plunge 



For the first grid point encountered on any blade 
other than the first, the computational molecule of figure 6 
is used. The new blade is assumed to protrude an 



inf initisimal 



distance past grid point F^ which is labeled 



F and F in the figure. Then, the equations developed 
22U 22L 

for grid points on the upper and lower surface of a blade 



are used to compute the perturbation quantities at F and 

. 22U 
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F 

22 L 



respectively. 



y 




x 



Figure 6. Computational molecule used for 
the first grid point on a new blade. 



(7) Wake Grid Poi nts 




Figure 7. Computational molecule used for 
a grid point in the wake. 



It is assumed that the first wake grid point is an 
inf initismal distance aft of the blade. To calculate tin 



48 



flow field quantities in the wake, the computational 
molecule is divided into two triangular molecules as shown. 
Then, the equations developed for a grid point on the upper 

surface of a blade apply to F while those for the lower 

22U 

surface of a blade apply to The com P uta tional 

molecule is shown in figure 7. 
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(III-59a) 

(III-59b) 

(III-59c) 

(III-59d) 



whe^e. 



K = U + C + A U 

12RU 1 1 R U 1 1 RU I 11IU 



K = U + C -AU 

i2iu mu mu i i iru 



(III-60a) 

(III-60b) 
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Eqs. (58) - (60) form a system of eight equations 
with twelve unknowns. However, in the linearized problem, 
the wake cannot support a pressure difference, and the 
normal velocity across the wake must be a constant. As a 
result of these two conditions, to first order effect, there 
will be no discontinuities reflected from the wake. 
Therefore, the relations, 

V = V = V (III-6 la) 
22RU 22RL 22R 



V = V = V 

22IU 22IL 221 



(111-6 1b) 



C 



22RU 



C 



22RL 



C 



22R 



(III-6 1c) 



C 



22IU 



C = C 

22IL 221 



(III-6 1 d) 



are used to reduce Eqs (58) and (59) to a determinant eight 
by eight system of equations which are: 
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(III-62a) 



u 


+ 


C 


_ 


A 


U 


K 




22RU 




22R 




I 


22IU 


12RU 




U 


+ 


C 


+ 


A 


U 


K 




22IU 




221 




I 


22RU 


12IU 




U 


+ 


V 


+ 


B 


< c 


u ) 


= K 


22RU 




22R 




I 


221 


22IU 


56R 


U 


+ 


V 




B 


( c 


u ) 


i = K 


22RU 




221 




I 


2 2R 


22RU 


561 



and 
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Egs. (60) remain unchanged. 
C. THE PRESSURE COEFFICIENT 



(III-62b) 
(III-62c) 
(III- 62d) 

(III-63a) 
(III-63b) 
(III-63C) 
(III-63 d) 



The non-dimensional 
surface of a blade is 
pressure difference of 
velocity perturbation 
Rewriting Eg. (11-48) as 

2 c 



pressure distribution along the 
determined using the linearized 
Eg. (11-48) and Teipel's sonic 
amplitude function. Eg. (21c) . 

- c 



(III- 64) 



and 



then dividing by the freestreara dynamic pressure yields 
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(Ill- 65) 




ikt 



Hence, the complex pressure difference across a blade at 
grid point F resulting from oscillation in pitch or plunge 



This expression for the pressure difference across a 
blade appears to be independent of interblade phase angle 8. 
However, the interblade phase angle was used for the 
computation of the normal velocity through the flow tangency 
equations, and they were used in the determination of C. 

Therefore, to compute the pressure distribution along the 
th 

n blade when it is in a specific position requires that 
the position be specified in terms of the position of the 
first blade. In general, the complex pressure is 



22 



is. 



P(x,y,t) = 



ikt 



?XpM 2 




(III-66a) 



ikt 



P (x,y ,t) = (P + iP ) e 
R I 



(III- 67) 



th 




blade is 



(kt) = (kt) + (n-1 ) 8 
n 



(III- 68) 



Then expanding Eq. (67) , 



P(x,y,t) = [P cos (kt) - P sin(kt) ] 



R 



I 



i[P sin(kt) + P cos(kt) ] 
R I 



(III- 69) 
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Figure 8 shows a complete oscillation cycle for a blade 
undergoing simple hamonic motion in pitch. The extension to 
plunge is obvious, and the following discussion is 
applicable there also. 




Figure 8. Blade oscillation cycle. 

The two conditions of primary interest here are when the 
th 

n blade is in the initial position 

(kt) = 0 (III- 70) 

n 

or in the mean position, pitching leading edge up, 

(kt) = 270 (III- 71) 

n 

Then solving for the position of the first blade when the 
th 

n blade is in the initial position gives, 

(kt) = - (n-1 ) 8 (III- 72) 

and in the mean position, 

(kt) = 270 - (n-1)S (III- 73) 

Except for flutter analysis, the real component of the 
pressure is normally of primary interest, and can be found 
by solving Eg. (69) with (kt) expressed as a function of 
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(kt) . For ease of showing the relation of the real 
n 

component of pressure for the two blade positions above, 
consider the special case of 8 = 0. Then, 

Re {P (x,y,t) } = P (III-74a) 

kt=0 R 

Re {P (x,y,t) } = P (III-74b) 

kt=270 I 

and thus. 

Re {P (x,y, t) } = Im (P (x,y,t) } (III-75a) 

kt=0 kt=270 



Re {P (x, y , t) } = Im {P (x,y,t) } (III-75b) 

1 ' s kt=270 1 kt=0 



Having obtained the pressure distribution over the 
blades, the dimensionless lift per unit span and the moment 
per unit span square are found from. 



L(x 



M (x 



,y,t) = = fp(x,y, 
J o 

.y.t) = = [( 

->o 



t) dx 



x-x ) P (x,y , t) dx 



(III-'76a) 



(III-76a) 



where L(x,y,t) and M(x,y,t) are complex quantities which 
include the effects of blade oscillation in pitch and 
plunge. Lift is positive upward and moments are positive 
clockwise (see figure 2) . 
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IV. FLUTTER ANALYSIS 



A. EQUATIONS OF MOTION FOR A SINGLE BLADE 

Lane [Ref. 14] showed that cascade flutter could be 

analyzed by considering only a single oscillating blade. 

Following Fung [Ref. 8], the equations of motion for a flat 

plate airfoil immersed in an inviscid flow are derived from 

a force balance of the inertia, elastic and aerodynamic 

forces. For bending motion, h is measured positive down from 

the elastic axis. The torsional motion is positive clockwise 

about the elastic axis (refer to figure 2) . It is assumed 

that a pair of springs, with spring constants K and K , 

a h 

oppose the blade motion. 



The total inertia force per unit span is, 

F = - (M h + S a) (IV- 1) 

i b b 

and this force exerts a moment about the elastic axis x of, 

o 

M = - (I a + S h) (IV- 2) 

i a b 

where 

M - total blade mass per unit span 
b 

S - wing static moment (about x ) 
b o 

I - wing mass moment of inertia (about x ) 
a o 

The springs produce restoring forces and moments of 

F = - hK (IV-3a) 

r h 
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M = " aK 



a 



(IV-3b) 



r 



Then the equations of motion for the blade are. 

Lift/span = M h + S a + hK 
b b h 



(IV-4a) 




(IV-4b) 



The assumption of simple harmonic motion is expressed as. 



where h and a are complex amplitudes. The blades 
o o 

oscillate with a constant interblade phase angle 8 . A phase 

difference exists between twist and flexure also, but cannot 
be determined until a flutter solution is obtained. Then, 
the ratio of blade deflection amplitudes will give this 
phase difference. 

With Eqs. (5) , the equations of motion for the blade 
become. 

Lift/span = - w 2 M h - W 2 S a + h X (IV-6a) 

b o b 0 0 h 

Moment/span = - o» 2 S h - W 2 I a + a K (IV-6b) 

b o a o o a 

The natural frequencies of twist and flexure are found by 
assuming that each mode can exist separately, in the absence 
of exciting forces. Then, 



i<nT 



h = h e 
o 



(IV-5a) 



icuT 



a = a e 
o 



(IV-5b) 




(IV-7a) 
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( IV— 7 b) 



tu = /K /I 
a a a 



K = Mu 2 
h h 



(IV-8a) 



and therefore. 



( IV-8b) 



lcuT 

Lift/span = [-h (a> 2 M - M <u 2 ) - a w 2 S ]e (IV-9a) 

o b b h 0 b 



Moment/span - [ -h cu 2 S - a (a> 2 I - I w 2 ) ]e (IV-9b) 

ob 0 a a a 



B. TWO DEGREE OF FREEDOM FLUTTER EQUATIONS 

The method used to investigate flutter in this report 
follows the procedure and notation used by Garrick and 
Rubinow [Ref. 9], who derived the lift force and moment 
equations, 

io)T h 

L = -4/>b 3 <u 2 e [ — (L + iL ) + a (L + iL ) ] 

G/R b 1 2 o 3 4 

(I V-10a) 

ieuT h 

M = -4flb 4 cu 2 e [ — - (M + iM ) + a (M + iM ) ] 

G/R H L b 1 2 o 3 4 J 

(IV- 10b) 

where L was defined as positive down. Equating the 

G/R 

magnitudes of the lift from Egs. (9a) and (10a) , one 
obtains. 
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—2 (L + iL - a + ft X) + a (L + iL - ax ) = 0 

b 1 1 2 ^ h ' 0 3 4 * a 



(IV-1 1) 



and by equating the amplitude of the moments. 



— (M + iM - ax ) + a (M + iM - ar 2 + SI X) = 0 
b 1 2 a 03 4 r a h 



(IV- 12) 



where the non-dimensional terms are defined as. 



(L + iL ) , (M + iM ) - real and imaginary component 

12 12 

of lift and moment due to plunge. 



(L + iL ) , (M + iM ) 



real and imaginary component 



4 3 4 

of lift and moment due to pitch. 

M. 

- wing density parameter ( 



4 pb* 

ratio of natural frequencies of flexure 



W. 2 

and torsion 



the quantity 

location of the wing center of gravity measured 
s b 

from the elastic axis ( — — ) 



M h b 



radius of gyration of the wing (- 



-) 



ratio of natural frequency in pitch to the 
t o 2 

circular frequency (— 3) 



Experimental results indicate the 
turbomachines £Ref. 18], 



in modern 



Eqs. (11) and (12) form an eigenvalue problem in ( — 

b 
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and (a Q ) which has a non-trivial solution if and only if the 

determinant of the coefficient matrix vanishes. Equating 
the determinant to zero yields a complex function in the 
unknown X, each component of which must be identically 
equal to zero. These equations are, 

ft ft X 2 + [ft (L - fi) + ft (M - ft ) ]X + C =0 (I V- 1 3a) 
ha al h3a fi 

(ft L + ft M ) X + C =0 (IV- 13b) 

a 2 h 4 I 

C = u.f x (M + L ) - (M -ft) - Lr 2 - llx 2 ] + D 
fi ^ a 1 3 3 a 1 a ^ a J fi 

(IV- 1 3c) 

C = tt[ x (M + L ) - M - L r 2 ] + D (IV-13d) 

I ^ a 2 4 4 2 a J I 

D = L M - L M - LM + L M (IV-13e) 

fi 1 3 3 1 2 4 4 2 



D = L h - L H + L M - L M (IV-13f) 

I 1 4 4 1 2 3 3 2 

For a fixed cascade 
geometry and interblade 
phase angle, Eqs. (13) 
are all transcendental 
functions of the reduced 
frequency k, or (1/k) , 
more accurately. The 
solution of Eqs. (13a) 
and (13b) is found using 
the Theodorsen method. 

/fie (X} and /I’m {X} are 
tabulated and plotted versus (1/k) . The flutter point is 
defined by the intersection of the real and imaginary curves 
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of /x” . The approximate shape of these curves is shown in 
figure 9. Since (1/k) appears transcendentally in every 
term of Eqs. (13)/ the curve of /Re {X} is not exactly 
parabolic, and for /im {X} , not actually linear. 



At the flutter point, the non-dimensional flutter 
frequency is. 




(IV- 14a) 



and the non-dimensional flutter speed is. 



Up _ 1 1 

oj c k/X 



(IV- 14b) 



k is defined by Eq. (11-39) . 



C. SINGLE DEGREE OF FREEDOM FLUTTER 



* For single degree of freedom oscillations in pitch, the 
eguation of motion for the blade is. 

Moment/span = I a + aK (IV- 15) 

a a 

which with the assumption of simple harmonic motion becomes, 

icuT 

Moment/span = [- w 2 I a + K a ]e (IV- 16) 

a o a 0 



The mechanical damping of the system is introduced by 

assuming the spring opposing the displacement in pitch to 

have a spring constant of the form K (1 + ig ) [Ref. 8]. 

a a 

The natural frequency in pitch is determined as before. 
Then, equating the amplitude of this expression for the 
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h 

moment with that of Eg. (10a) with (— ) = 0, one obtains, 

b 

ft ( X - 1 ) + M =0 (I V- 17a) 

a 3 

M + g ft X = 0 (IV- 17b) 

4 a a 

The solution method is the same as for the two degree of 
freedom problem, and the non-dimensional flutter frequency 
and speed are found from Eqs. (14). A similar development 
is required to obtain the equations for single degree of 
freedom plunge. However, based on the results of this, and 
other works, plunge oscillations are always stable. 

For both single degree of freedom and coupled 

th 

flexure- torsion motion, the flutter calculations for the n 
blade must be repeated for varying interblade phase angles 
until the minimum flutter speed is found. This speed is 
termed the critical flutter speed. For all speeds below it, 
cascade operation will be stable. 

An alternate method for determining the onset of flutter 
is to compute the real component (the in phase component) of 
the work per cycle of the blade on the freestream, 
r*%j 

Work/cycle = Re { / (M ctdt + L hdt)} (IV- 18) 

Jo “ h 

As long as the blade does positive work on the freestream, 
the cascade will be stable. However, if the aerodynamic 
damping decreases to the point that the airstream does work 
on the blade, flutter will occur. For single degree of 
freedom pitch oscillations, with zero structural damping, 
the sign of the imaginary component of the moment determines 
whether the work/cycle is positive or ne _tive, and can 
therefore be used to indicate the onset of flutter 
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